%% Series 6 - Exercise 2 for the course of DOE
%% Full factorial 2^4
% Table of Contents
% Full factorial
% Load data
% Analysis of the results
%     Check the matrix of experiments
%         Construction the term matrix (16 coefficients)
%         Labels of the 15 coefficients
%     Construction of the model matrix
%     Compute the half-effects
%         Bar chart of relative effects
%         Production of the normalplot
%         Half normal plot
%     Stepwise fit
%
%
% A process with 4 factors is tested with a full factorial design 2^4.
% The response is an output of the process whose sensitivity has to be determined.
% Compute the half-effecs and perform an ANOVA
%% Load data

data=readtable('data_doe_2020.xlsx','sheet','Series 4',"Range",'B3:G19');
disp(data)

%
%% Analysis of the results
%
%% Check the matrix of experiments
% The visual check of the matrix of experiments shows that it is a full factorial design (columns are constituted by +1 and -1 following fractions 1/2, 1/4, ...
%
%% Construction the term matrix (16 coefficients)
% The keyword "interactions" of the function fitlm() gives a linear model with only the interaction of first order. To have a full interaction model it is necessary to fed the function with a term matrix: one columns per factor and one row per coefficient.

modelspec=[0 0 0 0
    1 0 0 0
    0 1 0 0
    0 0 1 0
    0 0 0 1
    1 1 0 0
    1 0 1 0
    1 0 0 1
    0 1 1 0
    0 1 0 1
    0 0 1 1
    1 1 1 0
    1 1 0 1
    1 0 1 1
    0 1 1 1
  1 1 1 1];

Ncoef=size(modelspec,1)-1; % without the constant

%
%% Labels of the 15 coefficients
% Labels to be used in graphics
%

labels={'a_1';'a_2';'a_3';'a_4';...
    'a_{12}';'a_{13}';'a_{14}';'a_{23}';'a_{24}';'a_{34}';...
    'a_{123}';'a_{124}';'a_{134}';'a_{234}';...
    'a_{1234}'};

%
%% Construction of the model matrix
% The function x2fx() is used to compute the model matrix from the essay matrix.
%

X = x2fx(data{:,2:5},modelspec);

%
%% Compute the half-effects
% Computation of the half-effects with the least square algorithm
%

coefficients = X'* data{:,6} /16; %simplified algorithm LSF
const=coefficients(1);
effects=coefficients(2:16);

%
%% Bar chart of relative effects

figure()
bar(effects/const,'r')
set(gca, 'XTickLabel',labels)
title('Half relative effects','FontWeight','bold','FontSize',16)
xlabel('Coefficients')
ylabel('a_i/a_o')

%
%% Production of the normalplot
%
% The function normalplot() draws the function between the coefficients and the probability of the cummulative function. Putting a handle  (h=graphical_function) gives the possibility to change the properties of the objects within the graphic.
%

figure()
h=normplot(effects);
title('Normal plot','FontWeight','bold','FontSize',16)
xlabel('Sorted coefficients')
ylabel('norminv(p)- theoretical quantiles')

set(h(1),"MarkerSize",12,"Marker","o","MarkerFaceColor","red")
[E,index]=sort(effects);
text(h(1).XData+1, h(1).YData+.05, labels(index),'FontSize',16)

%
% Analysing the normalplot allows us to pretend that only a_4, a_{14} and a_{34} are significant.
% With a few coefficients less it is possible to perform an ANOVA :

modelspec2=[0 0 0 0
    1 0 0 0
    1 0 1 0
    0 0 0 1
   0 0 1 1
    1 0 0 1];

mdl1=fitlm(data{:,2:5},data{:,6},modelspec2)

%
%% Half normal plot
% Half normal plor is sometime easier to analyse. It can be produced with the routine probplot('halfnormal',x)

figure()
h=probplot('halfnormal',abs(effects));
title('Half normal plot','FontWeight','bold','FontSize',16)
xlabel('Sorted coefficients')
ylabel('norminv(p)- theoretical quantiles')

 [E,index]=sort(effects);
 text(h(1).XData+0.7, h(1).YData, labels(index),'FontSize',16)
 set(h(1),"MarkerSize",12,"Marker","o","MarkerFaceColor","red")

%% Stepwise fit
% An alternative is the stepwise :

mdl2=stepwiselm(data{:,2:5},data{:,6},modelspec)

% The stepwise method keeps more coefficients.
%
% The same result would have been obtained by limiting the model to the interactions of first order and eliminating coefficients with hight p-value.

mdl3=fitlm(data{:,2:5},data{:,6},'interactions')

modelspec4=[0 0 0 0
    1 0 0 0
    0 0 0 1
    1 0 1 0
    1 0 0 1
    0 0 1 1];

mdl3=fitlm(data{:,2:5},data{:,6},modelspec4)

%
%
% The final model is:
% y=11.9 + 5.9 x_1 - 14 x_4 -5 x_1x_3 +17.9 x_1x_4 -13.1 x_3x_4
% or eventualy adding an interaction of second order, what is correct statistically but probably without significant impact at the level of engineering $y=11.9 + 5.9 x_1 - 14 x_4 -5 x_1x_3 +17.9 x_1x_4 -13.1 x_3x_4 -3.9 x_1x_3x_4  $